Population class extension: distribution_functions module
Module containing the predefined distribution functions
The user can use any of these distribution functions to generate probability distributions for sampling populations
To add custom functions you can take any function and add it to the class instance before running the code. See https://stackoverflow.com/a/28060251 for some tips on how to do that
- There are distributions for the following parameters:
mass
period
mass ratio
binary fraction
- Tasks:
- TODO: make some things globally present? rob does this in his module..i guess it saves
calculations but not sure if I’m gonna do that now
TODO: add eccentricity distribution: thermal, Mathieu eccentricity
TODO: Add SFH distributions depending on redshift
TODO: Add metallicity distributions depending on redshift
TODO: Add initial rotational velocity distributions
TODO: make an n-part power law that’s general enough to fix the three part and the 4 part
- class binarycpython.utils.population_extensions.distribution_functions.distribution_functions(**kwargs)[source]
Bases:
object
Extension for the Population grid object that contains the distribution functions
- Arenou2010_binary_fraction(m)[source]
Arenou 2010 function for the binary fraction as f(M1)
GAIA-C2-SP-OPM-FA-054 www.rssd.esa.int/doc_fetch.php?id=2969346
- Parameters
m (
Union
[int
,float
]) – mass to evaluate the distribution at- Return type
Union
[int
,float
]- Returns
binary fraction at m
- Izzard2012_period_distribution(P, M1, log10Pmin=- 1.0)[source]
period distribution which interpolates between Duquennoy and Mayor 1991 at low mass (G/K spectral type <~1.15Msun) and Sana et al 2012 at high mass (O spectral type >~16.3Msun)
This gives dN/dlogP, i.e. DM/Raghavan’s Gaussian in log10P at low mass and Sana’s power law (as a function of logP) at high mass
TODO: fix this function
- Parameters
P (
Union
[int
,float
]) – periodM1 (
Union
[int
,float
]) – Primary star masslog10Pmin (
Union
[int
,float
]) – minimum period in base log10 (optional)
- Return type
Union
[int
,float
]- Returns
‘probability’ of interpolated distribution function at P and M1
- Kroupa2001(m, newopts=None)[source]
Probability distribution function for Kroupa 2001 IMF.
The (default) values for this is:
default = { "m0": 0.1, "m1": 0.5, "m2": 1, "mmax": 100, "p1": -1.3, "p2": -2.3, "p3": -2.3 }
- Parameters
m (
Union
[int
,float
]) – mass to evaluate the distribution atnewopts (
Optional
[dict
]) – optional dict to override the default values.
- Return type
Union
[int
,float
]- Returns
‘probability’ of distribution function evaluated at m
- Moe_di_Stefano_2017_multiplicity_fractions(options, verbosity=0)[source]
Function that creates a list of probability fractions and normalises and merges them according to the users choice.
TODO: make an extrapolation functionality in this. log10(1.6e1) is quite low.
The default result that is returned when sampling the mass outside of the mass range is now the last known value
Returns a list of multiplicity fractions for a given input of mass
- Moe_di_Stefano_2017_pdf(options, verbosity=0)[source]
Moe & diStefano function to calculate the probability density.
takes a dictionary as input (in options) with options:
M1, M2, M3, M4 => masses (Msun) [M1 required, rest optional] P, P2, P3 => periods (days) [number: none=binary, 2=triple, 3=quadruple] ecc, ecc2, ecc3 => eccentricities [numbering as for P above]
mmin => minimum allowed stellar mass (default 0.07) mmax => maximum allowed stellar mass (default 80.0)
- build_q_table(options, m, p, verbosity=0)[source]
Build an interpolation table for q, given a mass and orbital period.
m and p are labels which determine which system(s) to look up from Moe’s data:
m can be M1, M2, M3, M4, or if set M1+M2 etc. p can be P, P2, P3
The actual values are in $opts:
mass is in \(opts->{m} period is \)opts->{p}
Since the information from the table for Moe and di Stefano 2017 is independent of any choice we make, we need to take into account that for example our choice of minimum mass leads to a minimum q_min that is not the same as in the table We should ignore those parts of the table and renormalise. If we are below the lowest value of qmin in the table we need to extrapolate the data
Anyway, the goal of this function is to provide some extrapolated values for q when we should sample outside of the boundaries TODO: fix description to be correct for python
- calc_P_integral(options, min_logP, max_logP, integrals_string, interpolator_name, mass_string, verbosity=0)[source]
Function to calculate the P integral
We need to renormalise this because min_per > 0, and not all periods should be included
- calc_e_integral(options, integrals_string, interpolator_name, mass_string, period_string, verbosity=0)[source]
Function to calculate the e integral
We need to renormalise this because min_per > 0, and not all periods should be included
TODO: create function ot actually do this in a more general way
- calculate_constants_three_part_powerlaw(m0, m1, m2, m_max, p1, p2, p3)[source]
Function to calculate the constants for a three-part power law
TODO: use the power law_constant function to calculate all these values
- Parameters
m0 (
Union
[int
,float
]) – lower bound massm1 (
Union
[int
,float
]) – second boundary, between the first slope and the second slopem2 (
Union
[int
,float
]) – third boundary, between the second slope and the third slopem_max (
Union
[int
,float
]) – upper bound massp1 (
Union
[int
,float
]) – first slopep2 (
Union
[int
,float
]) – second slopep3 (
Union
[int
,float
]) – third slope
- Return type
Union
[int
,float
]- Returns
array of normalisation constants
- const_distribution(min_bound, max_bound, val=None)[source]
a constant distribution function between min=min_bound and max=max_bound.
- Parameters
min_bound (
Union
[int
,float
]) – lower bound of the rangemax_bound (
Union
[int
,float
]) – upper bound of the range
- Returns
returns 0
- Return type
returns the value of 1/(max_bound-min_bound). If val is provided, it will check whether min_bound < val <= max_bound. if not
- cosmic_SFH_madau_dickinson2014(z)[source]
Cosmic star formation history distribution from Madau & Dickonson 2014 (https://arxiv.org/pdf/1403.0007.pdf)
- Parameters
z – redshift
- Returns
Cosmic star formation rate in Solar mass year^-1 mega parsec^-3
- duquennoy1991(logper)[source]
Period distribution from Duquennoy + Mayor 1991. Evaluated the function self.gaussian(logper, 4.8, 2.3, -2, 12)
- Parameters
logper (
Union
[int
,float
]) – logarithm of period to evaluate the distribution at- Return type
Union
[int
,float
]- Returns
‘probability’ at self.gaussian(logper, 4.8, 2.3, -2, 12)
- fill_data(sample_values, data_dict)[source]
Function that returns the normalised array of values for given logmass and logperiod used for the e and q values
TODO: make sure we do the correct thing with the dstep
- flat()[source]
Dummy distribution function that returns 1
- Returns
1
- Return type
a flat uniform distribution
- flatsections(x, opts)[source]
Function to generate flat distributions, possibly in multiple sections
- Parameters
x (
float
) – mass ratio valueopts (
dict
) – list containing the flat sections. Which are themselves dictionaries, with keys “max”: upper bound, “min”: lower bound and “height”: value
- Return type
Union
[float
,int
]- Returns
probability of that mass ratio.
- gaussian(x, mean, sigma, gmin, gmax)[source]
Gaussian distribution function. used for e.g. Duquennoy + Mayor 1991
- Parameters
x (
Union
[int
,float
]) – location at which to evaluate the distributionmean (
Union
[int
,float
]) – mean of the Gaussiansigma (
Union
[int
,float
]) – standard deviation of the Gaussiangmin (
Union
[int
,float
]) – lower bound of the range to calculate the probabilities ingmax (
Union
[int
,float
]) – upper bound of the range to calculate the probabilities in
- Return type
Union
[int
,float
]- Returns
‘probability’ of the Gaussian distribution between the boundaries, evaluated at x
- gaussian_func(x, mean, sigma)[source]
Function to evaluate a Gaussian at a given point, but this time without any boundaries.
- Parameters
x (
Union
[int
,float
]) – location at which to evaluate the distributionmean (
Union
[int
,float
]) – mean of the Gaussiansigma (
Union
[int
,float
]) – standard deviation of the Gaussian
- Return type
Union
[int
,float
]- Returns
value of the Gaussian at x
- gaussian_normalizing_const(mean, sigma, gmin, gmax)[source]
Function to calculate the normalisation constant for the Gaussian
- Parameters
mean (
Union
[int
,float
]) – mean of the Gaussiansigma (
Union
[int
,float
]) – standard deviation of the Gaussiangmin (
Union
[int
,float
]) – lower bound of the range to calculate the probabilities ingmax (
Union
[int
,float
]) – upper bound of the range to calculate the probabilities in
- Return type
Union
[int
,float
]- Returns
normalisation constant for the Gaussian distribution(mean, sigma) between gmin and gmax
- get_integration_constant_q(q_interpolator, tmp_table, qdata, verbosity=0)[source]
Function to integrate the q interpolator and return the integration constant
- get_max_multiplicity(multiplicity_array)[source]
Function to get the maximum multiplicity. Used in the M&S calculations
- imf_chabrier2003(m)[source]
Probability distribution function for IMF of Chabrier 2003 PASP 115:763-795
- Parameters
m (
Union
[int
,float
]) – mass to evaluate the distribution at- Return type
Union
[int
,float
]- Returns
‘probability’ of distribution function evaluated at m
- imf_scalo1986(m)[source]
Probability distribution function for Scalo 1986 IMF (defined up until 80Msol): self.three_part_powerlaw(m, 0.1, 1.0, 2.0, 80.0, -2.35, -2.35, -2.70)
- Parameters
m (
Union
[int
,float
]) – mass to evaluate the distribution at- Return type
Union
[int
,float
]- Returns
‘probability’ of distribution function evaluated at m
- imf_scalo1998(m)[source]
From Scalo 1998
Probability distribution function for Scalo 1998 IMF (defined up until 80Msol): self.three_part_powerlaw(m, 0.1, 1.0, 10.0, 80.0, -1.2, -2.7, -2.3)
- Parameters
m (
Union
[int
,float
]) – mass to evaluate the distribution at- Return type
Union
[int
,float
]- Returns
‘probability’ of distribution function evaluated at m
- imf_tinsley1980(m)[source]
Probability distribution function for Tinsley 1980 IMF (defined up until 80Msol): self.three_part_powerlaw(m, 0.1, 2.0, 10.0, 80.0, -2.0, -2.3, -3.3)
- Parameters
m (
Union
[int
,float
]) – mass to evaluate the distribution at- Return type
Union
[int
,float
]- Returns
‘probability’ of distribution function evaluated at m
- interpolate_in_mass_izzard2012(M, high, low)[source]
Function to interpolate in mass
TODO: fix this function. TODO: describe the args high: at M=16.3 low: at 1.15
- Parameters
M (
Union
[int
,float
]) – masshigh (
Union
[int
,float
]) –low (
Union
[int
,float
]) –
Returns:
- Return type
Union
[int
,float
]
- ktg93(m, newopts=None)[source]
Probability distribution function for KTG93 IMF, where the default values to the three_part_powerlaw are: default = {“m0”: 0.1, “m1”: 0.5, “m2”: 1, “mmax”: 80, “p1”: -1.3, “p2”: -2.2,”p3”: -2.7}
- Parameters
m (
Union
[int
,float
]) – mass to evaluate the distribution atnewopts (
Optional
[dict
]) – optional dict to override the default values.
- Return type
Union
[int
,float
]- Returns
‘probability’ of distribution function evaluated at m
- linear_extrapolation_q(qs, indices, qlimit, qdata, end_index, verbosity=0)[source]
Function to do the linear extrapolation for q.
- merge_multiplicities(result_array, max_multiplicity, verbosity=0)[source]
Function to fold the multiplicities higher than the max_multiplicity onto the max_multiplicity
- if max_multiplicity == 1:
All the multiplicities are folded onto multiplicity == 1. This will always total to 1
- if max_multiplicity == 2:
The multiplicity fractions of the triple and quadruples are folded onto that of the binary multiplicity fraction
- if max_multiplicity == 3:
The multiplicity fractions of the quadruples are folded onto that of the triples
- number(value)[source]
Dummy distribution function that returns the input
- Parameters
value (
Union
[int
,float
]) – the value that will be returned by this function.- Return type
Union
[int
,float
]- Returns
the value that was provided
- poisson(lambda_val, n, nmax=None)[source]
Function that calculates the Poisson value and normalises TODO: improve the description
- powerlaw(min_val, max_val, k, x)[source]
Single power law with index k at x from min to max
- Parameters
min_val (
Union
[int
,float
]) – lower bound of the power lawmax_val (
Union
[int
,float
]) – upper bound of the power lawk (
Union
[int
,float
]) – slope of the power lawx (
Union
[int
,float
]) – position at which we want to evaluate
- Return type
Union
[int
,float
]- Returns
probability at the given position(x)
- powerlaw_constant(min_val, max_val, k)[source]
Function that returns the constant to normalise a power law
TODO: what if k is -1?
- Parameters
min_val (
Union
[int
,float
]) – lower bound of the rangemax_val (
Union
[int
,float
]) – upper bound of the rangek (
Union
[int
,float
]) – power law slope
- Return type
Union
[int
,float
]- Returns
constant to normalise the given power law between the min_val and max_val range
- powerlaw_constant_nocache(min_val, max_val, k)[source]
Function that returns the constant to normalise a power law
TODO: what if k is -1?
- Parameters
min_val (
Union
[int
,float
]) – lower bound of the rangemax_val (
Union
[int
,float
]) – upper bound of the rangek (
Union
[int
,float
]) – power law slope
- Return type
Union
[int
,float
]- Returns
constant to normalise the given power law between the min_val and max_val range
- powerlaw_extrapolation_q(qdata, qs, indices)[source]
Function to do the power-law extrapolation at the lower end of the q range
- raghavan2010_binary_fraction(m)[source]
Fit to the Raghavan 2010 binary fraction as a function of spectral type (Fig 12). Valid for local stars (Z=Zsolar).
The spectral type is converted mass by use of the ZAMS effective temperatures from binary_c/BSE (at Z=0.02) and the new “long_spectral_type” function of binary_c (based on Jaschek+Jaschek’s Teff-spectral type table).
Rob then fitted the result
- Parameters
m (
Union
[int
,float
]) – mass to evaluate the distribution at- Return type
Union
[int
,float
]- Returns
binary fraction at m
- sana12(M1, M2, a, P, amin, amax, x0, x1, p)[source]
distribution of initial orbital periods as found by Sana et al. (2012) which is a flat distribution in ln(a) and ln(P) respectively for stars * less massive than 15Msun (no O-stars) * mass ratio q=M2/M1<0.1 * log(P)<0.15=x0 and log(P)>3.5=x1 and is be given by dp/dlogP ~ (logP)^p for all other binary configurations (default p=-0.55)
arguments are M1, M2, a, Period P, amin, amax, x0=log P0, x1=log P1, p
example args: 10, 5, sep(M1, M2, P), sep, ?, -2, 12, -0.55
TODO: Fix this function! Half of the input here can be taken out and calculated within the function itself.
- Parameters
M1 (
Union
[int
,float
]) – Mass of primaryM2 (
Union
[int
,float
]) – Mass of secondarya (
Union
[int
,float
]) – separation of binaryP (
Union
[int
,float
]) – period of binaryamin (
Union
[int
,float
]) – minimum separation of the distribution (lower bound of the range)amax (
Union
[int
,float
]) – maximum separation of the distribution (upper bound of the range)x0 (
Union
[int
,float
]) – log of minimum period of the distribution (lower bound of the range)x1 (
Union
[int
,float
]) – log of maximum period of the distribution (upper bound of the range)p (
Union
[int
,float
]) – slope of the distribution
- Return type
Union
[int
,float
]- Returns
‘probability’ of orbital period P given the other parameters
- three_part_powerlaw(m, m0, m1, m2, m_max, p1, p2, p3)[source]
Generalised three-part power law, usually used for mass distributions
- Parameters
m (
Union
[int
,float
]) – mass at which we want to evaluate the distribution.m0 (
Union
[int
,float
]) – lower bound massm1 (
Union
[int
,float
]) – second boundary, between the first slope and the second slopem2 (
Union
[int
,float
]) – third boundary, between the second slope and the third slopem_max (
Union
[int
,float
]) – upper bound massp1 (
Union
[int
,float
]) – first slopep2 (
Union
[int
,float
]) – second slopep3 (
Union
[int
,float
]) – third slope
- Return type
Union
[int
,float
]- Returns
‘probability’ at given mass m